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Abstract 

The effects of strong Coulomb correlations in dense three-dimensional electron-hole plasmas are studied 
by means of unbiased direct path integral Monte Carlo simulations. The formation and dissociation of bound 
states, such as excitons and bi-excitons is analyzed and the density-temperature region of their appearance 
is identified. At high density, the Mott transition to the fully ionized metallic state (electron-hole liquid) is 
detected. Particular attention is paid to the influence of the hole to electron mass ratio M on the properties 
of the plasma. Above a critical value of about M = 80 formation of a hole Coulomb crystal was recently 
verified [Phys. Rev. Lett. 95, 235006 (2005)] which is supported by additional results. Results are related to 
the excitonic phase diagram of intermediate valent Tm[Se,Te], where large values of M have been observed 
experimentally. 
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I. INTRODUCTION 



Strongly correlated Coulomb systems have been in the focus of recent investigations in many 
fields, including dense plasmas in space and laborator3*i'2iM^ electron-hole plasmas in semicon- 
ductors, charged particles confined in traps or storage rings (see e.g.^ for an overview). In these 
systems the Coulomb interaction energy U is often larger than the mean kinetic energy K, i.e. the 
coupling parameter r = \{U)\/ {K) > 1. recently Coulomb and Wigner crystallization, which 
may occur when F ~ 100 ^ 1, has attracted much attention. Coulomb crystals were observed in 
ultracold trapped ions^'^'^, in dusty plasmas and storage rings^iiSiii. 

There exist many strongly correlated Coulomb systems where quantum effects are important. 
Examples are dense astrophysical plasmas in the interior of giant planets or white dwarf stars^ 
as well as electron-hole plasmas in solids, few-particle electron or exciton clusters in mesoscopic 
quantum dots (see^*^ and references therein). The formation of Coulomb bound states such as 
atoms and molecules or excitons and bi-excitons, of Coulomb liquids and electron-hole droplets^ 
exemplify the large variety of correlation phenomena that exist in these systems. 

Recombination of electrons and positive charges to neutral bound complexes strongly reduces 
the Coulomb coupling and thus acts against formation of Coulomb crystals in two-component 
charged particle systems. Quite recently the conditions for the existence of Coulomb crystals in 
neutral plasmas containing (at least) two oppositely charged components has been studied^^, 
and it was found that the mass ratio M of positive and negative carriers plays a crucial role. There 
exists a threshold value of M of about 80 in three-dimensional plasmas. Such values are possible in 
semiconductors, which leads to the prediction of hole crystallization in semiconducting materials 
with a sufficiently large effective mass asymmetr y ^^i^^'^^i^^ . Such values are feasible e.g. in the 
intermediate valence Tm[Se,Te] system, which under pressure and at very low temperature even 
might show the phenomenon of excitonic Bose condensation^-. 

In the present paper we substantially extend the analysis of previous work^^ on two- 
component partially ionized Coulomb systems with mass ratios M varying between one and about 
one thousand, thus covering plasmas, ranging from positronium, over condensed matter systems 
(almost) to hydrogen. Thereby we focus on the fundamental aspects of Coulomb correlations in 
two-component plasmas in dependence on M. Special emphasis is placed on situations with M 
close to the critical value for hole crystallization. 

From a theoretical point of view, the complex processes of interest which involve strong 
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Coulomb forces, quantum and spin effects are difficult to treat within the framework of analytical 
approaches. Therefore, over the last decade there has been a high activity in the development of 
numerical techniques capable to tackle strongly correlated Coulomb systems (plasmas 
A technique which is particularly well suited to describe equilibrium properties of two-component 
plasmas in the strong coupling and degeneracy regime, is the path integral quantum Monte Carlo 
(PIMC) method. Remarkable progress has been obtained in applying these techniques to Fermi 
systemsi'2i22i22i^. Since PIMC simulations of macroscopic Coulomb systems are hampered by the 
notorious fermion sign problem, several strategies have been developed to overcome or at least 
"weaken" this difficult y^^'^^^^^ . Within the restricted PIMC approach additional assumptions on 
the density operator were adopted, which reduce the sum over permutations to even (positive) 
contributions only^^^. This requires the knowledge of the nodes of the density matrix, however, 
which for interacting macroscopic systems are known only approximatel}*^. Hence the accuracy 
of the results, in particular in the regime of strong correlations, is difficult to assess. An alter- 
native approach are direct PIMC simulations which have occasionally been attempted by various 
groups, e.g.— but in general were not sufficiently precise and efficient for practical purposes. In 
recent years an improved path integral representation of the A^-particle density operator has been 
developed^^*^^ that allows for direct fermionic path integral Monte Carlo (DPIMC) simulations 
of dense plasmas in a large range of temperatures and densities, see"^^ for an introduction. 

The present paper applies the DPIMC method to electron-hole plasmas with strong mass- 
asymmetry. Here we consider situations where also the heavy component (referred to as "holes" 
hereafter) has to be treated quantum mechanically, unlike for hydrogen-like plasmas. A second ex- 
tension of our previous simulations is an improved treatment of the exchange-effects which allows 
us to reach higher densities and lower temperatures, as required to study the Mott effect and hole 
crystallization. Sec. |Il] gives a brief overview on the DPIMC approach for calculating thermody- 
namic quantities. Details on the derivations of the basic formulas, including equation of state and 
energy and a discussion of the quantum pair potentials used in the simulations are given in Ap- 
pendices |A] and IB An extensive numerical study of strongly correlated two-component Coulomb 
systems is presented in Sec. Hill Here we first give an overview on possible correlation effects in 
the limits of small and very large mass ratios. Then we consider more in detail the semiconductor 
system which is closest to the critical value of M = 80 for hole crystallization. In particular, 
energy and pressure, the microscopic electron-hole configurations, the fraction of bound states, 
as well as various (spin-dependent) pair distribution functions and charge structure factors were 
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calculated in a wide density and temperature range. Finally, numerical results for the hole crystal 
are presented. The main results are summarized in Sec. UVl 



II. PATH INTEGRAL MONTE CARLO PROCEDURE 

We consider a neutral two-component plasma consisting of Ne = = N electrons and holes 
in equilibrium with the Hamiltonian, H = K + If^, containing kinetic energy K and Coulomb 
interaction energy Lf^ parts. The thermodynamic properties at given temperature T and volume V 
are then completely described by the canonical density operator, p = e~^^ /Z, with the partition 
function 

Z{Ne,N,,V;P) = j^^Yl j dqpiq,^;P), (1) 

''■''■ay 

where P = 1/ ksT, and p(g, a; (3) denotes the diagonal matrix elements of the density operator at a 
given value cr of total spin 2;— projection. In Eq. ([T]), q = {ge, qh} and a = {o-g, ct/,} are the spatial 
coordinates and spin degrees of freedom of the electrons and holes, i.e. qa = {qi,a ■ ■ ■ (li,a ■ ■ ■ QNa,a} 
and (Ta = {(Ti a . . . a; a . . . (Tat^ „}, with a = e,p. All thermodynamic functions can be directly com- 
puted from the partition function. The resulting equations for density matrix, pressure (equation 
of state) and internal energy will be derived in Appendices lAl and iBl 

Of course, the exact density matrix of interacting quantum systems is not known (particularly 
for low temperatures and high densities), but it can be constructed using a path integral approach^ 
based on the operator identity e'^^^ = e"^'^^ ■ e~^^^ . . . e~^^^ , where A/3 = [3/{n + l). This 
allows us to express the density operator in terms of a product of (n + 1) known high-temperature 
density operators (at {n + l)-times higher temperature). In the coordinate representation this 
yields products of off-diagonal high-temperature density matrices (g*^'^^^^|e~^^^|g*^'^-') where k = 
1, . . . , + 1). Accordingly each particle is represented by a set of [n + 1) coordinates ("beads"), 
i.e. the whole configuration of the particles is represented by a 3(Ai'e + Nh){n + 1) -dimensional 

vecor ^ {,™, . . . ,<»i . . . ,<:;■', . c . . . ,;;«'}. 

Figure. [H illustrates the representation of one (light) electron and one (heavy) hole. The circle 
around the electron beads symbolizes the region that mainly contributes to the partition function 
path integral. The size of this region is of the order of the thermal electron wavelength \e{T), 
while typical distances between electron beads are of the order of the electron wavelength taken 
at an [n + l)-times higher temperature. The same representation is valid for each hole but it 
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electron 




FIG. 1: (Color onUne) Beads representation of electrons and holes. Here Ag = 2Tih^ [3 /m^, ^ = 
27r/i^A/3/me, qf'l = qf'} + Aa.c '7^ e\ and a = a' . The holes have a similar beads representation, however 
A^ /i is (m/i/me)-times smaller, so their beads distribution is not resolved in the figure. 



is not shown since, due to the larger hole mass, the characteristic length scales are substantially 
smaller. Nevertheless, in the simulations below, the holes are treated according to the full beads 
representation. Details, including the treatment of the spin, are given in the Appendices. 

To evaluate the density matrix, accurate results for the high-temperature approximation are 
necessary. As we have shown earlier^^, for sufficiently high temperature, i.e. for large number 
n of time slices) each high-temperature factor can be expressed in terms of two-particle density 
matrices (p = 1, . . . , A^^^, t = 1, . . . , A";,, a,b = e, h) 



pab{qp,a, q'p^a^ qt,b, q't^b'i P) 



[mamb 



,3/2 



X exp 



exp 



{lt,b - qt,b) 



exp[-/9<l> 



ab\ 



(2) 



with the familiar Kelbg potential^^ 



^ab{Xab](^) 



^ab-^ab 



^/nXab (1 - erf(Xafe)) 



(3) 



where Xn 



-7= C dte * . The derivation of 



'■ \Qp,a — (lt,b\/^ab, and the error function erf(x) 
approximation ^ and discussion of its accuracy is given in App.|Bl 

In our DPIMC scheme we use different types of steps, where either electron (hole) coordinates 
Qt,e (Qp.h) or individual electronic (hole) beads are moved until convergence of the calculated values 
is reached. Using periodic boundary conditions (PBC) the basic MC cell (filled yellow square in 
Fig. [21) is periodically repeated in x, y and z directions. 
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As mentioned above the main contribution to the path integral representation of the partition 
function comes from configurations for which the typical size of the clouds of electronic (hole) 
beads is of the order of the thermal wavelength of electrons (holes). At the moment our computer 
resources allow to consider up to about one hundred electrons and holes with several tens beads 
in the basic MC cell. Due to this limitation on the number of particles in the MC cell there is a 
restriction on the size of the MC cell for a given density. In the case of a highly degenerate plasma 




FIG. 2: (Color online) Sketch of the boundary conditions used for the simulations. Electron and hole clouds 
of beads are denoted by solid green and red lines respectively. The dashed line shows the cell for choosing 
the image of particle b which is closest to particle a. For further explanation see text. 

the thermal wavelength and the typical size of the electronic clouds of beads may be larger than 
the size of the basic MC cell. So beads of electrons belonging to the basic MC cell can penetrate 
into neighboring images of the main cell, as well as electronic beads from a neighboring cell can 
extend into the basic MC cell. Due to this fact we have improved our PIMC scheme compared 
to previous simulations: In our new calculations we assume that an electron (hole) belongs to a 
certain cell if its physical coordinate = 1, • • • , Na, a = e,h) belongs to this cell. In Fig. [21 
electron and hole clouds of beads are marked by solid green and red lines. Only a few periodic 
images of these particles are shown by related dashed lines. 

Let us consider the PBC for the calculation of pressure and energy in more detail. For distances 
between beads with the same number / in the Kelbg potentials and its derivatives [respectively first 
and second line in curly brackets in Eqs. (IB3I ) and (|B7I) 1 we used the standard PBC (see Fig. [2] a). 
Namely, in the Kelbg potential and its derivatives, instead of the distance ab between beads with the 
same number / of electrons i and j we take the smallest distance ab' to one of the electron images 
j'. The same applies to electron-hole and hole-hole distances. Furthermore, in calculations of the 
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scalar products and derivatives of the Kelbg potential [terms C and D in Eq. (IB3h and terms A and 
B in Eq. (IB7I) 1 the situation is more complicated due to the dependence of the scalar products on 
the angle between vectors to beads of the particle from the basic MC cell and its periodic images. 
In our calculations we first of all choose, for a given particle i, the nearest image j", according 
to the distance between coordinates and g°„ ^ only, as shown in Fig.[2]b). ab" is the smallest 
distance for all ab. Then for this pair i, j" we calculate all scalar product terms A, B, C, D and the 
related derivatives of the Kelbg potential. The same is done for electron-hole and hole-hole pairs. 

In our previous calculations determinants of the exchange matrices [cf. Eq. (IB6I) of App. |Bl 
were only computed for particles belonging to the basic Monte Carlo cell. However, with increas- 
ing degeneracy (n\^) the ratio of the particle thermal wavelength to the size of the Monte Carlo 
cell also increases. If this ratio approaches one exchange effects between particles in the main 
MC cell and their images in the neighbor cells have to be included. Therefore, in the present cal- 
culations we take into account the exchange interactions of electrons and holes from neighboring 
Monte Carlo cells, namely first from the (3'^ — 1) nearest-neighbor cells, then from the 5'^ — 1 
next nearest-neighbor cells and so on. These improved calculations were tested for both an ideal 
plasma and a non-ideal hydrogen plasma. Excellent agreement with the known analytical results 
for an ideal plasma was found up to densities where the parameter reaches values of several 
hundreds. 

In the present simulations of dense electron-hole plasmas, we varied both the particle number 
and the number of beads. We found that in order to obtain convergent results for the thermo- 
dynamic properties in the density-temperature range considered below it is sufficient to simulate 
systems with particle numbers of Ne = Nh = 50 . . . 100. Of course, the accuracy is strongly 
affected by the number of beads n. To exclude an n-dependence of our calculations, the density 
matrices in the high-temperature decomposition were always taken at temperatures above the exci- 
ton binding energy. In practice, a number of about n = 20 beads turn out to be sufficient. In order 
to simplify the computations further, we included only the dominant contribution in the sums over 
the total electron and hole spin, which corresponds to s = Ne/2 electrons and k = Nh/2 holes 
having spin up and down, respectively. The contribution of the other terms is small and vanishes 
in the thermodynamic limit. Let us emphasize that for all results presented below the maximum 
statistical error is about 5%, which is sufficient for the present analysis. Note that this accuracy 
can be achieved at an acceptable cost of computer time. Of course, the error can be systematically 
reduced by increasing the length of the Monte Carlo run. 
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III. SIMULATION RESULTS 



We now apply the theoretical scheme developed in the preceding sections to a partially ionized 
dense electron-hole plasma. We will be interested in strong Coulomb correlation effects such as 
bound state (excitons, bi-excitons, clusters), their modification by the surrounding plasma and 
their eventual breakup at high densities due to pressure ionization (Mott effect). Beyond the Mott 
density, we expect the possibility of hole crystallization if the hole mass is sufficiently large^. 
To detect these effects, we have extended our first-principle DPIMC simulations to a large range 
of mass ratios, temperatures and densities. Below, the density of the two-component plasma is 
characterized by the Brueckner parameter, = d/as, defined as the ratio of the mean distance 
between particles d = [3/47r(ne + rih)]^/'^ and the exciton Bohr radius = h^e/e^nij., where rif. 
and Uh are electron and hole densities, respectively, and e is the background dielectric constant. 
In what follows we compute and discuss the spatial particle configurations, the pair distribution 
functions, static structure factors, fractions of electrons and holes in bound states, the internal 
energy and the equation of state. 

A. Electron hole plasma with small/large mass ratio M 

Let us first discuss Coulomb correlations in two-component plasmas for the limiting cases of 
small and very large mass ratios. In Figs. l3l4l we show data for M = 1 (positronium, left columns 
in the two figures) and M = 952 (right columns, illustrating the situation typical for hydrogen and 
plasmas of other chemical species). In addition, we consider the cases of small and large density, 
corresponding to = 10 (upper row) and = 0.33, (lower row in Figs. [3] and U) respectively, 
and temperatures well below the the exciton binding energy Eb = /2eaB- To allow for a direct 
comparison with the Tm[Te,Se] system below, we fix the effective electron mass to me = 2.1mo 
(mo is the free electron mass) and the background dielectric constant to e = 25 which leads to a 
binding energy EB/kB = 517 K being two orders of magnitude smaller than for positronium and 
hydrogen. 

In Fig. [3] we present typical spin-resolved "snapshots" of the electron-hole state in the sim- 
ulation box for small and large M (left and right columns, respectively). One clearly sees the 
influence of the hole mass on the particle probability distribution (corresponding to the DeBroglie 
wavelength): In the left part, electron and hole "clouds" have the same size, whereas in the right 
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FIG. 3: (Color) Snapshots of the electron-hole plasma configuration in the simulation box for two mass 
ratios, left (right) column: M = l,T = WOK (M = 952, T=50K) and two densities corresponding to, 
upper (lower) row: = 10 (r^ = 0.33). Spin up and spin down electrons (holes) are marked by yellow 
and blue (red and pink) clouds of dots, respectively, and the Monte Carlo cell is given by the gray grid lines 
(PBC were used). 

part, holes have practically shrunk to dots (negligible size compared to the interparticle distance). 
Let us first concentrate on the low-density limit, = 10, upper row in Fig. |3l At this density 
Coulomb correlations are strong enough to give rise to bound states - excitons. They are clearly 
visible in the snapshots from pairwise clustering of electrons and hole "clouds " for both, small and 
large mass ratios. Occasionally, also clusters of three or four particles are seen which correspond 
to trions (exciton ions) and bi-excitons, respectively. 

A quantitative analysis of the behavior is obtained by computing the pair distribution functions 
(PDF) which are shown in the upper row of Fig. |4l They are computed in the simulations from the 
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density operator according to 

9ab{r) = / dqS{ri^a- qi,a)S{r2M- q2,b) p{q,cr]p) . (4) 
" V 

For both mass ratios we observe a clear peak of r'^geh at about r^h = Iob, corresponding to 
excitons. For the large mass ratio, Fig. |4b), the peak is substantially higher which is a consequence 
of the increase of the binding energy by a factor two compared to the case M = 1 (the reduced 
mass increases from me/2 to mg) which stabilizes the bound states. This trend is also seen in the 
snapshots (upper panel): The fraction of electron "clouds" closely attached to holes is significantly 
higher in the right figure. Further, the h-h PDF for the large mass ratio. Fig. Sb), show signatures 
of bi-exciton formation with a peak distance of about lAas- Also the e-e PDF shows peaks, one 
at lower distances, correspponding to electrons with different spin projections located between the 
holes and a weaker pronounced one at larger distances. These peaks are not seen for M =1. 

Let us now turn to the limit of high densities, lower row in Figs. [3] and IH Here the influence 
of the mass ratio is even more dramatic, leading to a qualitative change of the plasma behavior. 
While the electrons are practically delocalized over the whole simulation volume for both M, 
the behavior of the holes changes from delocalized (small mass ratio; lower left part of Fig. [3]) 
to fully localized (large mass ratio; lower right part). Obviously, in both cases no bound states 
exist. Instead, we observe a Fermi gas-like state of electrons and holes, at small Af , and a hole 
crystal which is embedded into an electron Fermi gas, at large M. This interpretation is confirmed 
by the behavior of the PDF (lower row of Fig. H]). They are almost structureless at small M. In 
contrast, at large M, pronounced peaks at finite r are visible in the hole-hole PDF being typical 
for a crystalline structure^. 

Thus, Figs. |3] and |4] allow us to conclude that in order to observe a crystal of charged particles 
(hole crystal) in a two-component quantum plasma, three requirements have to be fulfilled: (i) 
sufficiently low temperature, (i) sufficiently high density (which causes pressure ionization of the 
bound states) and (iii) a large mass ratio. Below, these requirements will be studied more in detail. 

The most interesting question is the role of the mass ratio. Crystals of ions (e.g. nuclei of 
carbon and oxygen) in the presence of a quantum electron gas are commonly accepted to exist in 
astrophysical objects such as White Dwarf stars'^, where the mass ratio is of the order M ~ 10^. 
But ion crystallization is expected to be possible also for much smaller mass ratios: proton crys- 
tallization (M = 1836) in dense hydrogen has been found in our previous simulations^^, see also 
Ref. and we have also found crystals of alpha partilces in pure heliuroi^. Figs. [3] and |4] indicate 
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FIG. 4: (Color online) Pair distribution functions (PDF) for the four combinations of mass ratio and density 
shown in Fig.[3l i.e., Fig. a): M = 1, = 10, b): M = 952, = 10, c): M = 1, = 0.33, d): M = 952, 
r<j = 0.33. Full black lines: e-e PDF, dashed red lines: h-h PDF, dash-dotted blue lines: e-h-PDF, dotted 
black line: r'^f^gehi'f'eh)- Distances are in units of the exciton Bohr radius as, note the different scales in the 
upper and lower row. 

that crystallization of holes should be possible even for M below 1000. Two-component plasmas 
with M < 1000 exist in condensed matter systems, such as semiconductors. In most traditional 
semiconductors, however, typical values of M are 3 ... 20. There have been predictions of the 
possibility of M ~ 100 in CuCl or Bi-Sb alloys under pressure, e.g.-^^, but without experimental 
evidence so far. 



B. Thermodynamic properties of electron hole plasmas with M = 40 

The largest mass ratios in condensed matter systems were experimentally reported by Wachter 
and co-workers^^ for the intermediate valent Tm[SexTei_x] alloys under pressure. In these mate- 
rials f-d-hybridization provides a narrow dispersive f- valence band and, as a consequence, a large 
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effective hole mass of the order of 50-100 (bare) electron masses. This system is of particular 
interest because of the long life time of the electron-hole plasma. Moreover the mass ratio is close 
to the predicted critical value of M ^ 80 for hole crystallization^-. 

TmSeo.45Teo.55 at ambient conditions is an indirect semiconductor with a gap of = 
130 meV. An excitonic level has been observed with Eb ~ 50 — 70 meV below the bottom of 
the d-band. Applying pressure the gap can be tuned (and even closed), and the material is spec- 
ulated to realize in the pressure region between 5 and 1 1 kbar an excitonic insulator, at least at 
very low temperatures^, the search for which has been run for a long time3 (see the experimental 
phase diagram in Fig. [5]). A necessary pre-condition is the existence of a large number density of 
(up to 10^° — 10^^ per cc) excitons of intermediate size (in order to avoid too strong overlap of the 
excitonic bound states)^^. For pressures exceeding 11 kbar exciton ionization has been observed 
as a result of the Mott effect. This is the region in phase space where hole crystallization should 
occur. 

In the following, we analyze this interesting system more in detail performing direct PIMC 
simulations within the parabolic band approximation. While this is certainly a strongly simplified 
model, we expect that the main trends, related to density and temperature variation, will be repro- 
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FIG. 5: Excitonic phase diagram of TmSeo.45Teo.55 taken from Ref.— . Measured points are designated by 
symbols. "Isobars" in the semi-conducting and semi-metallic phase are marked as dotted lines, whereas the 
"isobar" entering the excitonic phase is denoted by a full line. The lower abscissa gives the corresponding 
energy gap AE (here negative values refer to the metallic state). The left boundary (at high pressure) of the 
exciton-rich phase corresponds to the pressure ionization (Mott effect). For more details see Ref.-^^^. 
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duced. The simulations were performed for an e-h plasma with rUe = 2.1mo ,'mh = 80mo, and 
e = 25. Several characteristic temperatures within the interval of the measurements by Wachter et 
al.— (T < 300i^) are chosen, which are well below the exciton binding energy. 

1. Pressure and internal energy 

We first consider the behavior of pressure p and internal energy E which are determined in 
our simulations exploiting formulas (IB7I) and (IB3I) . derived in the App. El Figure [6] shows p 
and E versus the Brueckner parameter (i.e. versus density to the power —1/3) at several fixed 
temperatures. At low densities and high temperature (300K), p and E reflect the classical ideal 
gas behavior. Reduction of density or/and temperature lead to a significant deviation from this 
behavior: Now p and E decrease due to the (overall attractive) Coulomb interaction in the plasma 
and, in particular, due to bound state formation of excitons and biexcitons, reaching a minimum 
around ~ 2 ... 4. For higher densities, p and E start to increase again monotonically. This 
is triggered by quantum effects - the plasma behaves as a Fermi gas. At the same time, bound 
states are expected to break up as a consequence of the Mott effect. The behavior of bound states 
will be verified below from the snapshots and pair distribution functions. The negative values of p 
observed at the lowest temperature point towards an instability of the homogeneous plasma state 
against formation of droplets or clusters. Electron-hole droplet formation in semiconductors is 
well established and was observed experimentally three decades ago^. This effect is similar to 
the so-called plasma phase transition discussed by many authors for dense hydrogen and other 
plasmas (see^i^^i^^iiME and references therein). 

2. Particle configurations 

Figure |7] shows snapshots of the electron-hole configuration in the simulation box at different 
temperatures and densities. According to the temperature decomposition of the density matrix, 
each electron and hole is represented by several beads. We used n = 20 beads, so for the upper 
panels, at T = 50 K, the high-temperature density matrices are taken at a temperature of lOOOii', 
being two times larger than Eb- The spatial distribution of the beads of each quantum particle is 
proportional to its spatial probability distribution. Figure Vindicates that for M 40 the typical 
size of the cloud of beads for electrons is several times larger than the one for the holes. Note, that 
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FIG. 6: Isotherms of pressure and internal energy versus Brueckner parameter (r<j ~ n 

-1/3) for an e-ri- 

plasma with m-g = 2.1mo, rrih = 80mo and e = 25. For comparison, also results for e = 10 at T = 300K 
are shown. Pressure is in units of the classical ideal pressure, P^f = {ue + nh)kBT. 



in the present strongly correlated system, the extension of electron and hole probability densities 
maybe quite different from the deBroglie wavelength which corresponds to the ideal case. 

At low density (r^ = 10, left column) and temperature T = 50/^, practically all holes are 
closely covered by electron beads, which means that electrons and holes form bound states, and the 
e-h plasma consists mainly of excitons. This interpretation will again be supported by the behavior 
of the pair distribution functions discussed below. Raising the temperature at fixed density leads 
to a (temperature-induced) ionization of the bound states. As a result we find a substantial number 
of free electrons and holes in the simulations (the degree of ionization will be shown in Fig. [TTl 
below). 

For intermediate densities (r^ = 4, middle column) we observe the formation of bi-excitons and 
many-particle clusters at low temperatures (T = 50K). Now the electron-hole system is strongly 
inhomogeneous. In this case the mean distance between particles d is of the order of the electron 
wavelength. At T = 200K bi-excitons and many-particle clusters are absent due to thermal break 
up of exciton-exciton bound states. 

At high density (r^ = 1) the electron wavelength exceeds the mean inter-particle distance d and 
even the size of the Monte Carlo cell used in the simulations, which is seen by the large extension 
of the electron beads clouds. At the same time excitons become unstable because two electrons 
bound to neighboring holes start to overlap allowing for electron tunneling from one exciton to 
another (pressure ionization, Mott effect), and the system transforms into a plasma of electrons and 
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FIG. 7: (Color) Snapshots of particle configurations at three particle densities corresponding to = 10 
(left), r<j = 4 (center), and = 1 (right). Spin up and spin down electrons (holes) are marked by yellow and 
blue (red and pink) clouds of dots, and the Monte Carlo cell is given by the gray grid Unes (PBC were used). 
Upper and lower panels show results at low (T = 50K) and high (T = 200K) temperatures, respectively. 

holes. Since the hole wavelength is significantly smaller than the electron wavelength (and may 
still be smaller than rf), in a certain region of -values the structure of the hole beads resembles a 
liquid state (lower right panel). 

3. Pair distribution Junctions and structure factors 

Now let us discuss the behavior of the spin-averaged electron-electron-, hole-hole- and 
electron-hole PDFs defined in Eq. dH). Figure [8] shows the three functions gab{r) versus the inter- 
particle distance r = ^ — r2,b for the same densities and temperatures as in Fig. |71 Here no 
smoothening has been carried out, i.e. the fluctuations of gee, Qhh and Qeh at small r reflect the 
magnitude of the statistical errors of our simulation. 

As an effect of the Coulomb and Fermi (statistics) repulsions, gee and ghh are suppressed at 
small distances. Due to the large mass difference the decay of the e-e correlations, however, is 
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essentially different from that of the h-h correlations. The asymptotic of Qee at small distances is 
determined mainly by electrons with opposite spin projections (since the strong Fermi repulsion 
is absent). For these electron pairs the main contribution to the repulsion comes from the effective 
quantum Kelbg potential (recall that it is finite at zero distance), allowing for tunneling of electrons 
up to zero separation. This is supported by the behavior of the spin dependent pair distribution 
shown in Fig. [lOl Of course, the tunnelling effect are much more pronounced for the lighter 
electrons. 

The strong peak of Qeh at low densities (r^ = 10 . . . 8) is caused by excitons. This is confirmed 
by considering the function r'^gehir) which exhibits a pronounced maximum at about las, i.e. at 
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FIG. 8: (Color online) Pair distribution functions ^ee (black solid line), g^h (red dashed line), and g^h (blue 
dot-dashed line) at T = 50K (left two columns) and T = 200-fir (right two columns). The dotted lines 
show r'^geh where, for better visibility of the exciton peak, the data are divided by a factor 120 in panels a 
and b, by a factor of 30 in figures c, d, g, and h, and by a factor of 10 in panels e, f, and i-1. 
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the exciton Bohr radius. At these densities, the functions gee and ghh exhibit no peak structure, 
i.e. there is no indication of formation of bound exciton-exciton complexes (such as bi-excitons or 
electron-hole droplets). With increasing temperature more and more excitons dissociate, and the 
maximum of Pgeh is reduced. 

Varying, at low temperatures, the density over three orders of magnitude, the maximum of r'^geh 
is suppressed and finally vanishes. At about Vg = 6, there is evidence that recombination into bi- 
excitons takes place. At the same time, the position of the maximum in r'^geh shifts from 1 to 2- 
3 ttB approximately, indicating the increasing radii of the bound states. This scenario is confirmed 
by the behavior of gee and ghh' For intermediate densities (2 < < 6), they exhibit distinct 
peaks at larger r, pointing towards the formation of bi-excitons and many-particle clusters. For 
Tj, < 1, the fraction of excitons is further reduced due to many-body effects (pressure ionization), 
and bi-excitons and many-particle clusters vanish. 

Finally let us relate the width of the peak of geh to the extension of the ground state wavefunc- 
tion (~ las)- At low densities (r^ > 6) the peak of geh is rather broad, indicating the population 
of excited states, while at high densities (r^ ~ 1) the width of the lowest peak in geh becomes 
significantly smaller than the corresponding width of the ground-state peak. At the same time the 
hole-hole pair distribution function reveals an ordering of the holes into a fluid-like state. 

Figure [9] shows the variation of the static charge structure factor, defined in reciprocal fc-space 

as 

^^^^^ ^ drr^iQabjr) - 1) sin(fcr)/ (fcr) 
\ drr^{gab{r) ~ l)\ 

According to Eq. ^ positive (negative) values of Sabik) indicate attraction (repulsion) in momen- 
tum space, see also Ref.— . 

Starting at low densities (r^ = 10), the negative values of the e-e and h-h structure factors at 
small momenta (large distances) originate from the strong Coulomb repulsion of quasi-free equally 
charged particles. The maximum of Seh, as well as the modulations of See and Shh at finite k, are 
indicative of exciton formation. Excitons set a new length scale in the structure factor at about 
0.2 [1/20^]. As expected, these signatures are washed out at higher temperatures, where excitons 
break up. At high density (r^ = 1) the large mass ratio between electrons and holes is responsible 
for the different behavior of gee and ghh- This is especially true for low temperatures when the 
electrons are perfectly delocalized but the holes still have fluid-like short-range correlations as 
noted above. The value of the normalization constant (denominator) in Eq. ^ is responsible for 
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FIG. 9: (Color online) Static structure factor See (black solid line), Shh (red dashed line), and Seh (blue 
dot-dashed Une) at T = 50K (left columns) and T = 200K (right columns). 

the different magnitude of the Shh at T = 50K and T = 200K. 

The e-e- and h-h PDF and structure factors presented in Figs. [8] and |9] were averaged over the 
spin degree of freedom of the particles. Our DPIMC simulations, however, allow direct inspection 
of the spin effects as well. Especially for the case of large degeneracy (n\^ ^ 1) we expect 
qualitatively different behavior of gH and glj^ at small distances, due to the influence of the Fermi 
statistics (a = e, h). This is confirmed by Fig. \W\ In other words, we clearly see the exchange- 
hole known from the ideal Fermi gas, but here it appears in a strongly interacting e-h plasma. In 
contrast, for the much heavier holes the quantum statistical repulsion is much less pronounced, 
and the decay of ghh for r — > is mainly triggered by the (spin-independent) Coulomb interaction. 

The snapshots depicted in Fig.|7]have shown that in the density regime 2 < < 6 the plasma 
contains a large fraction of bi-excitons. Now, Fig. [TOl confirms this conclusion (consider, e.g., the 
case = 4): There is not much difference in (7^^ and gjj^ - both functions exhibit a similar peak 
at a hole-hole distance of about Sa^. In contrast, the electron behavior is strongly spin dependent: 
There is a high probability to encounter two electrons with different spin projections between two 
holes (at a distance la^ . . . 20^) - which is the expected behavior known from bi-excitons (or the 
hydrogen molecule). 

The corresponding charge structure factors clearly show that the difference between electrons 
and holes with respect to the spin correlations are most pronounced in the regime where bound 
states are formed but, of course, the differences are less important at very small and large k (corre- 
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FIG. 10: (Color online) Spin resolved pair distribution functions gah{f) and structure factor Sah{k) for 
electrons (solid black lines) and holes (dashed red lines) at T = 50 K. Left (right) columns correspond to 
antiparallel (parallel) spin projections of the particle pair. 

sponding to large and small distances, respectively). We found attraction (repulsion) for opposite 
(equal) spin projections for both particles as A; ^ 0, except for the case of high densities, where the 
hole fluid forms. Here, in the vicinity of the Mott transition, the difference between the structure 
factor for electrons and holes with same and opposite spin projections is notably smaller; they are 
close to their spin-averaged values given in Fig. |9]for T = 50 K and = 2. 



4. Degree of ionization. Mott density 

The above analysis of the snapshots and PDFs has indicated that, in a broad range of den- 
sities, the e-h-plasma is partially ionized containing a varying fraction of free and bound elec- 
trons and holes. The value of the critical Mott density (corresponding to the Brueckner parameter 
^Mott-j \Yhere excitons vanish due to pressure ionization is essential for the occurrence of hole 
crystallization. In particular, the critical mass ratio was found to directly depend on r^°** as 
j\,fcr _(_ 1 = r^'^/r^"**, where r^e = and rf ^ 100 is the critical value for Coulomb crys- 
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tallization in a one-component plasma . The value of r^°** depends in a complicated way on the 
structure of the excitons (which may be quite different for different materials) and on many-body 
effects such as screening and bound-state renormalization. This is a very complex problem which 
has been discussed in a variety of approximations the accuracy of which, however, is difficult to 
assess (see e.g. Ref.^ for an overview). Based on our DPIMC simulations avoiding any simpli- 
fying assumptions, we have the possibility to obtain a consistent many-body result for r^°". To 
this end, we first need to find a definition for the degree of ionization which is applicable to our 
simulations. 

The main problem of PIMC simulations in configuration space is that there is, in principle, no 
clear subdivision into bound and free "components" possible. Electron-hole correlations arising 
from bound and scattering states contribute to the same quantities, such as energy, equation of 
state or PDF. Note that the usual subdivision in energy space into bound and scattering states 
according to negative and positive relative pair energies Ea also breaks down in the case of strong 
correlations: Eigenvalues Ea and wavefunctions ^Q,(r) are being renormalized, \E'a ^a', Ea — > 
Ea- In the vicinity of the Mott point, the distinction between bound and scattering states becomes 
meaningless, bound state levels merge into the scattering continuum''. 

Nevertheless, a rough estimation of the fraction 1 — a'"" of e-h bound states (a;*°" denotes the 
degree of ionization) can be obtained by analyzing the PDF geh- In particular, at low temperatures 
ksT <^ Eb and far below the Mott density, the plasma will consist of excitons in the ground state 
El, i.e. geh{r) ~ |^i(r) ^ (see^^). In the general case, the PDF will contain a superposition (mixed 
state) of all renormalized eigenfunctions, weighted with the Boltzmann factor. 



While the scattering states are delocalized, the bound state wavefunctions are localized in space 
leading to an increase of g^^h beyond the result for an ideal plasma (7*^ = 1. Due to the normaliza- 
tion of geh values of geh > 1 at small e-h distances must be compensated by a depletion, geh < 1, 
at larger distances. Now, recall that geh is related to the probability density by Pehir) ~ r'^geki'^) 
which, in the case of excitons, is strongly enhanced around reh = Io-b, cf. Fig.[8l Hence we can 
use as a "pragmatic" definition of the fraction of bound states (this idea is due to N. Bjerrum who 
used it very successfully in the theory of electrolytes^) the probability of e-h-pairs being at small 




(6) 



a 
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distances (where Qeh > 1): 




= 1 — a 



ion 



lo [9eh{r) - 1] dr 



(7) 



Jo r^geh{r)dr 



Here is the second zero of r'^{geh — 1) which is located to the right of the exciton peak. By 
subtracting dr in the nominator, the uncorrelated contributions are eliminated from the full 
probability density. The denominator accounts for the normalization (giving the full probability 
of particles being in bound or scattering states). Of course, this definition is only qualitative since 
it does not exclude attractive scattering states and thus may slightly overestimate the "true" bound 
state fraction. 

An analogous ratio can be defined for the hole-hole bound-state (bi-exciton) fraction 



with [r*', Tfo] being the interval where (ghhir) — 1) is positive. 

We expect that the definition (|7]) is well suited for a numerical determination of the Mott density 
as the density, where the plasma becomes dominated by free particle behavior. We will identify 
j^Mott ^YiQ density where the bound state fraction falls below 5 . . . 10%, i.e. = 90 . . . 95%. 

Figure [m presents our DPIMC results for the e-h bound state fraction according to Eqs. 
([8]). At low densities (r^ > 6) and low temperatures (T < 100 K) practically all electrons and 
holes are bound in excitons, i.e. bi-excitons and many-particle clusters are absent. Increase of 
temperature, T > lOOi^', leads to a strong ionization of excitons, i.e. free particles dominate. 
In the intermediate density regime (2 < < 6) and for low temperatures (T < 100 K), the 
largest fraction (up to about 20%) of h-h bound states (bi-excitons and many-particle clusters) is 
observed. It is interesting to note that at higher temperatures (T ~ lOOK. . . 300 K) the exciton 
fraction increases a bit for higher density. 

A decrease of the bound state fraction below 5% is observed for ^ 1 which we, therefore, 
identify with the Mott density, r^°** ^ 1 with an error of about 30% (which corresponds to 
r^e ~ 1-2 which was used in Ref.'^). From this result, we confirm the value of the critical mass 
ratio (see above) M''^ ^ 83 which is in the range of previous predictions, e.g.^^. Since also the 
value for rf is expected to have an error of about 20%, a total uncertainty of about 50% has 
to be expected. This means that hole crystallization might occur already below M ^ 80 which 
underlines again the interest in the Tm[SexTei_x] system. 




hh 



jib' r'^ [ghh{r) - 1] dr 
XI' r'^ghhir) dr 



(8) 
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FIG. 1 1 : Fraction of the electron-hole (top left, a)) and hole-hole bound states (top right, b)) versus Brueck- 
ner parameter for different temperatures and dielectric constants. Bottom Fig. (c) illustrates the temper- 
ature induced break up of excitons for different mass ratios and dielectric constants. 

Let us, therefore, come back to the experimental phase diagram. Fig. [51 There the excitonic- 
rich phase should exist up to a temperature of about 250K. Then it is of course interesting to 
calculate the temperature dependence of the e-h bound state fraction in this regime (cf. the results 
for Tg = 10, bottom part Fig.fTTI). Although our simulations do not exhibit a sharp transition from 
an excitonic to an ionized "phase", we can again formally use a bound state fraction of 10 % as 
the boundary of the exciton dominated phase. Then, at low density, e.g. for = 10, these bound 
states will be stable up to T = 150 . . . 250 K. Note that these values are quite sensitive to e which 
is not precisely known. We therefore have performed simulations for two slightly different values. 
Increase of e leads to a more weakly bound excitons and consequently reduces the ionization 
temperature. The same tendency is observed if the mass ratio M is lowered. Fig. [TT] shows that 
the ratio of temperatures corresponding to a*"" = 50 % is the same as the corresponding ratio of 
the reduced masses = menih/ {rUe + rrih) (keeping e fixed). For example, a ratio of about 2 is 
observed for the data belonging to rrih = 80 and = 2.1, and 1 .5 for those belonging to rrih = 80 
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and ruh = 4.2. Clearly, as expected for an isolated exciton, the binding energy is proportional to 
rrir. Note however that in the present case the excitons are embedded into an e-h plasma which 
influences the bound state spectrum, so the recovery of the single exciton behavior is not a trivial 
result. 

We can also give a rough estimate of the critical electron density where the transition from 
an insulator (built up of e-h bound states) to a (semi-) metal takes place, again using a value of 
Q,«o7i _ ggo^ ^ criterion. The critical value of ~ 1 corresponds to a particle density of the 
order of 10^^ cm^'^ which is in reasonable agreement with the estimates from the experiments. 

C. Increased hole localization with increasing mass ratio 

After we have confirmed the pressure induced break up of excitons and determined the Mott 
parameter r^°** in Tm[SexTei_x] with a mass ratio of M ^ 40, we can now return to the question 
of strong h-h correlations at high density and the possibility of hole crystallization. To this end we 
have performed DPIMC simulations at low temperatures and high densities, well above the Mott 
point, and varied the mass ratio over two orders of magnitude. 

Fig. [T2l shows the results for Vs = 0.5 and ksT / Eb = 0.064 with M ranging from 5 to 400. To 
simplify the analysis we do not mark electrons and holes having different spin projections with dif- 
ferent colors (we found that, for hole crystallization, spin effects are of minor importance). Since 
the density is an order of magnitude higher than the Mott density, electrons and holes are in the 
plasma state. For all values of M the electrons are completely delocalized, forming a Fermi gas 
which is very weakly correlated. At the same time, the nature of the hole state changes drastically 
increasing of M. Initially (M = 5), the holes are also in a Fermi gas-like state where individual 
holes penetrate each other. An increase of M up to 50 leads to continuous growing hole localiza- 
tion. At M = 50 individual holes can be distinguished: Their DeBroglie wavelength becomes 
smaller than the average distance between two holes dh- A further increase of M to 100 leads 
to an additional reduction of Xh by a factor -\/2 (with dh unchanged) and, consequently, to hole 
localization. The holes form a regular lattice - a Coulomb crystal (bottom left Fig.). Increasing 
the mass ratio further, the lattice becomes more rigid (bottom right Fig.), until the holes become 
practically point-like (cf. Fig.©. Based on these simulation results we can conclude that, at these 
values of density and temperature, hole crystallization occurs between M = 50 and M = 100, 
which is in surprisingly good agreement with the analytical estimate for M^^ based on the Mott 
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FIG. 12: (Color) Snapshots of the electron-hole plasma at high density (r^ = 0.5) and T/ Eb = 0.064. The 
pictures show the spin-averaged electron (yellow) and hole (red) beads configuration for different values of 
the mass ratio M (from top left to bottom right): 5, 12, 25, 50, 100, 400. 

density obtained above for the same parameters. 

Obviously, the snapshots allow for a very rough determination of the critical mass ratio only. A 
more accurate estimate can be obtained by analyzing the relative distance fluctuation of the holes 
in dependence on M. These results, reported in Ref.^^, confirm the above result for M"' . We also 
have studied the thermal properties of the hole crystal which arise to a temperature dependence of 
^Mott gggj^ from Fig. [m (left part) the bound state fraction isotherms has reached saturation 
at about 60K (corresponding to slightly more than O.I-Eb)- Thus our result for M^'^ will not 
change significantly at lower temperatures, thus reflecting the ground state behavior. The full 
phase diagram of the hole crystal was published in Ref.— . 
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IV. DISCUSSION 



In this paper we have presented a numerical (computer simulation) analysis of strong Coulomb 
correlations in dense three-dimensional two-component plasmas at rather low temperatures. We 
were in particular interested in systems with a mass ratio M ~ 40 which is intermediate be- 
tween "usual" semiconductors (M ~ 2 ... 10) and "conventional" plasmas (such as hydrogen 
with M > 1836). Two-component Coulomb systems with this mass ratio behave in many aspects 
like plasmas, with the main difference that quantum and spin properties of the heavy particles 
(holes) cannot be neglected. Such systems might be realized in intermediate valence semiconduc- 
tors under pressure. 

From a plasma physics standpoint these are very interesting materials because they allow to in- 
vestigate nontrivial high density phenomena of current interest, such as pressure ionization (Mott 
effect), plasma phase transition and metal-insulator transition. These effects exist in "conven- 
tional" plasmas only above a density of the order of Ig cm^'^ which, meanwhile, can be achieved 
in laser or ion beam compression experiments, however, only for short periods of time which 
makes precision measurements of the plasma properties very difficult. The same physical effects 
can be observed in the above mentioned semiconductor materials under stationary (equilibrium) 
conditions at densities and temperatures which are easily accessible in experiments. In fact, the 
(qualitative) phase diagrams of dense plasmas and electron-hole plasmas are readily translated 
from one to another by rescaling the binding energy and the Bohr radius^^. For this reason, the 
phase diagram of Tm[Se,Te], shown in Fig. [51 is of interest also for the physics of dense partially 
ionized plasmas. 

In this paper, we concentrated on the central part of this phase diagram, where a large fraction 
of bound states exists, on the Mott transition and on the high density effects (hole crystallization). 
A theoretical treatment of these effects is very difficult, because many-body effects such as bound 
state formation, screening and quantum effects have to be taken into account self-consistently. 
While for such complex systems analytical methods fail, our first principle path integral Monte 
Carlo approach is well suited. In order to avoid additional approximations which are particularly 
questionable in the region of the Mott point, we have used direct fermionic simulations. With an 
improved treatment of the spin statistics we were able to present reliable simulations. Restricting 
the simulations to temperatures above 6% of the exciton binding energy and densities to values not 
significantly higher than the Mott density allowed us to avoid serious difficulties arising from the 
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fermion sign problem. 

Most notably, we have shown that, above the Mott point, two-component plasmas with large 
mass anisotropy show interesting Coulomb correlation phenomena: with increasing density holes 
can undergo a phase transition to a Coulomb liquid and to a Wigner crystal which are embedded 
into a degenerate electron Fermi gas. Such crystals are expected to exist in White Dwarf stars 
where the mass ratio exceeds 10^. However, crystal formation in a two-component plasma should 
be possible also for the light elements, such as hydrogen^^ '^'' and helium and, more generally, for 
plasmas with mass ratios as low as 80 (cf.'^). This should be possible to achieve in the semi- 
conductor materials announced in this paper. More subtle questions, such as the symmetry of the 
crystal and its energy cannot yet been answered conclusively, mainly because of the yet too small 
size of the simulations (50 electrons and 50 holes are presently feasible). Therefore, in order to 
obtain more accurate data, e.g. for the internal energy of a macroscopic two-component plasma at 
very high density, a significant increase of the simulation size would be highly desirable. 
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APPENDIX A: PATH INTEGRAL REPRESENTATION OF THERMODYNAMIC QUANTI- 
TIES 



Let us consider a neutral two-component plasma consisting of N^, = = N electrons and 
holes in equilibrium with the hamiltonian, H = K+lf^, containing kinetic energy K and Coulomb 
interaction energy U'^ = U^h ~^ ^ee + ^e/i contributions. The thermodynamic properties in the 
canonical ensemble with given temperature T and fixed volume V are fully described by the 
density operator p = e^^^ /Z with the partition function ([T]). 

Pressure and internal energy follow from 



13E 



dlnZ _ 
dV ~ 
_ d\nZ 



a dlnZ 

W da 



a=l 



dp 



(Al) 
(A2) 



where a = L/LqIS a length scaling parameter. 

The density matrix of interacting quantum systems is can be constructed using a path integral 
approach^ based on the operator identity e~^^ = e~^^^ ■ e~^^^ . . . e~^^^, where AjS = P/{n + 
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1), which allows us to rewrite the integral in Eq. ([T]) 

^ Pe Ph 

The spin gives rise to the spin part of the density matrix (S) with exchange effects accounted for 
by the permutation operators Pg and Ph acting on the electron and hole coordinates and spin 

projections a'. The sum is over all permutations with parity Kp^ and Kp,^. In Eq. (|A3I) the index k = 
1 . . .n + 1 labels the off-diagonal high-temperature density matrices p^'^^ = p [q^''^^\ g^'^^; A/?) = 
|/g(fc-i)|g-A/3//|q,(fc)^^ Accordingly each particle is represented by a set of + 1 coordinates 
("beads"), i.e. the whole configuration of the particles is represented by a 3{Nf, + Nh){n + 1)- 
dimensional vector q = {q[^l, . . . . . . . . . q^^^^^; qf^ . . . q^^^}^} (^^^ ^^S- 

To determine the energy in the path integral representation (|A3I) each high-temperature density 
matrix has to be differentiated in turn (here we extend our earlier hydrogen results of Refs.— i^): 



n+1 

X > p-'...p 

fc=l 



/3 



O- Pe Ph 

dp(k) 



dp 



(A4) 



)("+!) =o{0).cr' 



One can show that the matrix elements p'^^^ can be rewritten as 



^(fc-l)|g-A/3H|^(fc)^ 



|p{A.))(p{A.)|e-A/3i^|^ 



• • • k^'^) 



(A5) 



where p('=)(p(^)) are conjugate variables to q^^ ^\q^''^). Evaluating the derivatives in Eq. (|A4|) 
further, it is convenient to introduce dimensionless integration variables rj^^^ = i^e^), where 



(k) 



Ka{q^a^ — q^a ) for A; = 1 , . . . , 77,, and = nia/ (27r/i^A/?) = 1/A^ ^. The main advantage 



is that differentiation of the density matrix now affects only the interaction terms 



(A6) 



where 
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with pi*^^ = p^a^/{Kah), p{k) = {pl^\pi^^), and use has been made of Eq. (|A5I) . Furthermore 

X(0) ^ (/t.gf , X^^) ^ (xf \ Xi'^)) with xi'^ = K^qT + Eti ^i'^ runs from 1 to n), 

and X("+i) = {kJ^^^\ njr^^^) = X(o). 
For = n + 1, we have 



hPfS 



Then, together with Eq. (|A4I) . we obtain for the energy 

= ^(Ar^ + AT;,) - i— i 

e 



V 



{ n+l 



X 



I. fc=i 



(A8) 



(A9) 



This way the derivatives of the density matrix have been calculated and we turn to the next point: 
Finding approximations for the high-temperature density matrices p^'^^ . 



APPENDIX B: HIGH-TEMPERATURE ASYMPTOTICS FOR THE DENSITY MATRIX. 
KELBG POTENTIAL 

In this section we discuss approximations for the high-temperature density matrix that can be 
used for efficient DPIMC simulations. This involves effective quantum pair potentials $^5, which 
are approximated by the Kelbg potential (see also previous work—). 



1. Pair approximation and Kelbg potential 

The N-particle high-temperature density matrix is expressed in terms of two-particle den- 
sity matrices (higher order terms become negligible at sufficiently high temperature, i.e. for 
large number n of time slices) given by (O. This results from factorization into kinetic and 
interaction parts, pab ~ Po Pab^ which is exact in the classical case, i.e. at sufficiently high 
temperature. The error made at finite temperature vanishes with the number of time slices as 
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l/m? (cf. Ref.— ). The off-diagonal density matrix element Q involves an effective pair inter- 
action which is expressed approximately via its diagonal elements, $^/^(gp_a, g^^, qt^., q[ ^; (3) ^ 
[^ab{<lp,a — <lt,b] /3) + ^abil'p^a " 'it fei for wMch wc usc the familiar Kelbg potentiaP^-^^ ©. 

Note that the Kelbg potential is finite at zero distance which is a consequence of quantum effects. 
The validity of this potential as well as of the diagonal approximation is restricted to temperatures 
substantially higher than the exciton binding energy.-^Siil which puts another lower bound on the 
number of time slices n. For a discussion of other effective potentials, we refer to Refs.'^^'^^'^Siil. 

Summarizing the above approximations, we can conclude that with the approximations ([21), 
^ each of the high-temperature factors on the r.h.s. of Eq. (IA3I) . carries an error of the order 
l/(n + 1)^. Within these approximations, we obtain the result 

pW = pfe~^^^^^''-''^{X^^-^^ - XW) + 0[(l/n + \f\ , (Bl) 

(k) 

where Pq is the kinetic density matrix, and U denotes the sum of all interaction energies, 
each consisting of the respective sum of pair interactions given by Kelbg potentials, [/(X^'^)) = 



2. Estimator for the total energy 

Let us now return to the computation of thermodynamic functions and derive the final expres- 
sions which follow from Eq. (IBll) and which will be used in the simulations. First, we note that in 
Eq. (IA9I ), special care has to be taken in performing the derivatives of the Coulomb potentials with 

respect to {3: Products {3 — have a singularity at zero inter-particle distance which is 

integrable but leads to difficulties in the simulations. To assure efficient simulations we transform 
the e-e, h-h and e-h contributions in the following way: 

d 



da I rfx('=-i)(X('^-i)l^-^'^"^ 



(A/3f/^(X(^-i); 



X (xC^-Dje-^^^^-^^^lX'^) + O {l/{n + if) 



+ 0(l/(n+l)2) 



(B2) 



This means, within the standard error of our approximation, O (n ^), we have replaced the sum of 
the Coulomb potentials by the corresponding sum of Kelbg potentials U, which is much better 
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suited for MC simulations. This result coincides with expressions, which can be obtained if we 
first choose an approximation for the high-temperature density matrices p'^'^^ using Kelbg potential 
and then take the derivatives. 

Thus, our final result for the energy is 
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and '^ab{x 



Apd[P'(^ab{x-, P')]/dP'\i3i=/:^p. Further, (. . . | . . . ) denotes the scalar product, and qpt, Vpt and Xpt 
are differences of two coordinate vectors: Qpt = qp,h - qt,h, Tpt = qp,e - qt,e, Xpt = qp^e - qt,h. 
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The density matrices psk appearing in Eq. (IB3I) are given by 
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where U = ^—^ Yl {^^^(^i'^ + Uh{xi;\AP) + f/,,(xi^ Xf, A/5) } , (B5) 

^ 1=0 

and (f)\ = exp[— Trir/f -* p], 0^ = exp[— 7r|r/p p]. Notice that the density matrix (IB4I) does not contain 
an explicit sum over the permutations and thus no sum of terms with alternating sign. Instead, the 
whole exchange problem is contained in the following determinant which is a product of exchange 
matrices of electrons (index s) and holes (index k) where s (k) denotes the number of electrons 
and holes having the same spin projections (or more details, we refer to Ref.^-), 
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3. Estimator for the equation of state 



In similar way, expressions for all other thermodynamic functions can be derived. Here, we 
only provide the corresponding result for the equation of state, 

PpV ^ 1 (3Z)-i 



s=0 k=0 



X ■ 



p=it=i d\q,t\ ^ d\r,,\ 



ddetM 



a oaeT.\\'(ljpf\\sk 
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(B7) 



:h B{x,,) = A{rl,) = and = 

Eqs. (IB3I) for the total energy and (IB7I ) for the pressure are readily understood: the first terms 



on the r.h.s. correspond to the classical ideal gas part. The ideal quantum contribution, in excess 
of the classical one, plus all correlation contributions are contained in the integral terms. The 
Coulomb correlation contributions arise from the terms with the Kelbg potentials ^ab, whereas the 
exchange contributions arise from the derivatives of the exchange matrix (last term). 

While there exist many alternative representations of the thermodynamic functions, the main 
advantage of the present expressions Eqs. (IB3I ) and (IB7I ) for energy and pressure is that the explicit 
sum over permutations has been converted into the spin determinant which can be computed very 
efficiently using standard linear algebra methods. Furthermore, each of the sums in curly brackets 
in Eqs. (|B3I) and (IB7I) is bounded when the number of high-temperature factors increases (n 
oo). Note that Eqs. (IB3I) and (IB7I) contain the important limit of an ideal quantum plasma in a 
natural wav^^. 
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